///*
// * ReactionPath.cpp
// *
// *  Created on: 2 Aug 2011
// *      Author: Allan
// */
//
//#include "ReactionPath.h"
//
//// C++ includes
//#include <fstream>
//#include <iomanip>
//#include <set>
//
//// GeoReact includes
//#include "KineticSolver.h"
//#include "EquilibriumSolver/EquilibriumSolver.h"
//#include "Utilities.h"
//#include "Equilibrator.h"
//
//ReactionPath::ReactionPath() :
//T(25.0), P(1.0)
//{}
//
//ReactionPath::ReactionPath(const vector<string>& allSpecies, const vector<string>& eSpecies, const vector<string>& kSpecies, const vector<string>& initialConditions) :
//T(25.0), P(1.0)
//{
//	Initialiaze(allSpecies, eSpecies, kSpecies, initialConditions);
//}
//
//void ReactionPath::Initialiaze(const vector<string>& allSpecies, const vector<string>& eSpecies, const vector<string>& kSpecies, const vector<string>& initialConditions)
//{
//	multiphase = CreateMultiphase(allSpecies);
//	
//	eReactions = CreateReactionSystem(eSpecies, multiphase.GetSpeciesMap());
//	
//	kReactions = CreateReactionSystem(kSpecies, multiphase.GetSpeciesMap());
//	
//	this->initialConditions = initialConditions;
//	
//	initialEquilibriumConditions = ExtractEquilibriumConditions(initialConditions, multiphase, eReactions);
//}
//
//void ReactionPath::SetTemperaturePressure(double T, double P)
//{
//	this->T = T;
//	this->P = P;
//}
//
//void ReactionPath::Solve(double tInit, double tEnd, double tDelta, const string& fileName)
//{
//	n = GuessSpeciesComposition(initialConditions, multiphase, T, P);
//	
//	SolveInitialEquilibrium();
//	
//	// Print the initial state of the multiphase system
//	cout << "Initial State of the Multiphase System: " << endl << endl;
//	
//	multiphase.PrintState(T, P, n);
//	
//	KineticSolver kineticsSolver(multiphase, eReactions, kReactions);
//	
//	kineticsSolver.SetControlSolution(1.0E-6);
//	
//	ofstream file(fileName);
//	
//	file << scientific << setprecision(8);
//	
//	file << setw(18) << left << "Time";
//	file << setw(18) << left << "DeltaTime";
//	BOOST_FOREACH(const string& species, multiphase.GetSpecies()) file << setw(18) << left << species;
//	file << endl;
//	
//	double t = tInit, dt = tDelta;
//	
//	while(t < tEnd)
//	{
//		cout << "Progress: " << 100.0 * t/tEnd << endl;
//		
//		file << setw(18) << t;
//		file << setw(18) << dt;
//		for(int i = 0; i < n.size(); ++i) file << setw(18) << left << n[i];
//		file << endl;
//		
//		kineticsSolver.Iterate(t, tEnd, dt, T, P, n);
//	}
//}
//
//void ReactionPath::SolveInitialEquilibrium()
//{
//	// The number of primary species
//	const unsigned Nj = eReactions.Products().size();
//	
//	// Check if the number of initial equilibrium conditions equals the number of primary species
//	if(initialEquilibriumConditions.size() == Nj)
//	{
//		Equilibrator equilibrator(multiphase.GetSpecies(), eReactions.Reactants(), initialEquilibriumConditions);
//	
//		equilibrator.SetTemperaturePressure(T, P);
//	
//		equilibrator.Solve(T, P, n);
//	}
//	else
//		cout << "Warning: The number of initial equilibrium conditions is not enough to determine the initial equilibrium state. Proceeding anyway." << endl;
//}
